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TECHNICAL PAPER 


A SCHEME FOR BANDPASS FILTERING MAGNETOMETER MEASUREMENTS 
TO RECONSTRUCT TETHERED SATELLITE SKIPROPE MOTION 


I. INTRODUCTION 


The lirst Eight oi NASA s tethered satellite system is presently scheduled lor 1492. It is to be 
carried onboard the space shuttle to a 300-km orbit with a 28.5° inclination 1 1 1. There the satellite will be 
deployed outward along the orbit radius vector, to a tether length reaching 20 km. and later retrieved. It is 
well known that the tether may exhibit some very strange dynamics, one of which is the so-called 
skiprope motion |2|. This mode ol vibration is of particular concern because it has very little damping 
and could pose a problem during satellite retrieval and docking. If this happens, the present plan is to 
have the shuttle (light crew execute an orbiter yaw' maneuver that is properly phased relative to the rotat- 
ing tether so as to cause the rotary motion to decay quickly [3|. However, to know when to maneuver, at 
what rate, and in what direction, requires knowing tether position as a function of time: and this means 
estimating it. The baseline approach is to do this on the ground, in near real-time, using a state observer. 
Originally, the observer used only satellite rate gyro measurements as inputs [4], Afterwards, it was 
modified to use satellite magnetometer measurements also. 

This paper describes a completely different approach to this problem. It uses the same magneto- 
meter measurements, but filters them through bandpass filters tuned to the skiprope frequency. It is 
simple to implement, quite robust, and can be used in parallel with or as a backup to the observer. It is 
presented in this paper in the following manner. A theoretical development of its underlying equations is 
given in section II. These are utilized in the implementation scheme described in section III. Simulation 
results lor it are given in section IV. Final comments are made in section V. 


II. THEORETICAL DEVELOPMENT 


Previous studies indicate a relationship between tether position caused by skiprope motion and 
satellite attitude, as shown in figure 1 |5|. Tether deflection at its midpoint, d. is related to the list of the 
satellite i sy axis from the orbit radius vector. a sw , by the relationship 



TT 


(I) 


where L is tether length. As the tether^kipropes. the satellite i sy axis cones about the orbit radius vector. 
iw . at the skiprope frequency with the /<,•* axis lying in the plane of the tether, as shown in figure I . Hence, 
il the coning motion can be reconstructed, the tether motion can be interred. The scheme proposed here 
seeks to do that based on measurements of the Earth s magnetic field with an onboard magnetometer. 



i w : Orbit Radius Vector 



Figure 1. Tether skiprope motion. 


To develop this scheme, two coordinate frames are utilized. One is the satellite body-tixed trame 
(7 and the other is the orbit-fixed frame (q.iUiK shown in figure 2. As noted there. t< is 
aliened with the orbit velocity vector. K with the orbit north, and i w with the orbit radius vector directe 
awav from the Earth. The <£,.&.&> frame is referenced toAe OWv.M frame by the Euler angles 
(a ' a .., a . . > The Earth's magnetic field flux density vector B is referenced to it by the parameters . 
fr /.Pi/). as shown in figure 3. where B is its magnitude in Gauss. It is assumed that the steady-state 
attitude motion of the satellite due to tether skiprope can be described by 


ct S 2 = a .s:v/ *i n ( w .sV) 
a s\ = CX.SIA/ C0S • 


where u, w is the skiprope frequency and can be either positive or negative depending on the direction of 
coning. The angle a,., is assumed constant and small enough, along with a, : and a,,, that cos (a Sf ) - 1 
and sin <a v ,) = a v/ . for / = 1.2.3. 
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Given this model, it is straightforward to show that 
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and 


g = -fi COS P/.;/ sin P ,/ /r + B cos Pat/. cos P '/ *V + B s7 " P" ' l 


From equation 


ations (3) and (4). the components of the Earth's magnetic field in satellite axes 


B si = -B cos P/-7 sin P,z + «.s? B cos Paa cos p.«-(B sin Pa/.) «s: 
g v , = B cos Pa/, cos p. A/ + av? B cos Paa sin Paz+ ( B sin P f.i) a vi 
B sy = B sin Paa -(B cos Pax cos p A/ > a.si ~(B cos Pax sin P'^ 1 • 


By the same 
equation (2) 


token, these are the satellite magnetometer outputs resolved into satellite axe; 
into equation (5) and assuming for the remainder of the development that 


a S M = Ot.SI.W — a S2M • 


vields. after some manipulation. 


_ _ B s i + B cos Paa (sin P^~«s3 cos $_ A / ± ^ sin ^ cos {(J)SR t + Tt/2) 

Bs\c g~~ 

„ Bsz~B cos Paa (cos P.az + a sx s * n P A ^ = aci , s jn B r , cos (a> S 7 ?f) 

«s:c g 

= B.V3-B sin Pax = cos [ Ws *,- (p. AZ + ir)l 
sx B cos Paa 


(4) 

are given by 


/ <5> 
i. Substituting 


I (6) 
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Hence, if the magnetometer outputs in satellite axes. (fi S i .B s2 .B sy ). are modified according to equation 
(6). the results are. in theory, three harmonics from which the following inferences can be made about 
satellite attitude motion, and hence tether motion, due to the skiprope phenomenon: 

( 1 1 ot s u is determined by the amplitude of B syc . Tether deflection at the midpoint of the tether can 
then be determined using equation (1) and knowing the tether length L. 

(2) The period of B S \( . T SK . determines the magnitude of u> s7f by the relationship 


<- *>SK = - » T sh . 


If fl.sic leads B s:c , then w v/? > 0. If B slc lags B S2C . then oi SR < 0. 

^ (3) When B syc is_^a maximum, then m S r ’ = $.\? + v. At this point. T sy lies in the plane formed by 

iw and B. tilted toward B from i w by the angle ctsu- Then i sy at any other point in time is readily deter- 
mined. assuming [3 \/ und u> S y? are known. Figure 4 helps to clarify this. There. B P and i S \ P are the projec- 
tions of B and i sy . respectively, onto the /(• — /\ plane. Tether position is it rad out of phase with £ v >. 

Hence, equation (6) and these associated inferences form the theoretical basis for the scheme proposed in 
this paper. A practical implementation of it is presented in section III. 



Figure 4. Projection ol i s2 and B onto /(- — /\ plane. 
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III. IMPLEMENTATION SCHEME 


In implementing equation (6). it is assumed that the satellite magnetometer outputs are teleme- 
tered to the ground and resolved into satellite axes, if necessary, to generate (fi s i .B S 2 -B S %). The Earth s 
magnetic field parameters (fi.Bzv .Ba/> are estimated on the ground via a math model of the Earth's mag- 
netic field and knowledge of the satellite orbit. Since a S \ is small, the terms in it are discarded to obviate 
the need to estimate it^and thus simplify the scheme. Consequently, equation (6) is implemented as 


fit i fi cos (if/ sin B \y , • o , , . -) t 

B su = — — — = a U / sin B«. cos (ta SK t + it/-) 

fi 

fi s -.-fl cos B/ t COS B \y . o . 

B s , c = — 1 - a.v.v/ sin {3, / cos (u> v/f /) 

fi 


fiv \ ~ B sin B/ / . i , /a , _,i 

B S \t = — = ot.v.i/ cos |u>sr/ (paz + 7T M • 

fi cos Bf ,. 


(7) 


where C) denotes an estimated variable. If the only motion ot the satellite was that caused by tether 
skiprope. then equation (7) could be utilized directly. However, there will be satellite motion due to other 
effects that create unwanted signals in (Bs\,Bs 2 -Bsi) an< ^ hence (fi.sic-fi.s' 2 ofi.s' 3 C'b An example ot this is 
pendulous motion, which is satellite oscillation about its center of mass. This is expected to have ampli- 
tudes on the order of several degrees at a frequency in the neighborhood of 0.03125 Hz. Motion due to 
tether skiprope is expected to be on the order of degrees or fractions of a degree, depending on the tether 
length |5). Its frequency is around 0.005 Hz for tether lengths greater than 2 km. To eliminate unwanted 
signals like this. (Bs\c-Bs 2 cBsx) ' n equation (7) are each filtered by a bandpass filter whose transfer 
function has the form 


G(.v) 


24 is/*},,) 

(i'/u),,)" + 2£(s/u>„) + 1 


where £<< 1 and co„ = <o s * ideally. At oj = co„ = uj s/? . the filter passes the skiprope frequency with 
unity gain and zero phase shift while attenuating all other frequencies. Hence, it filters higher frequency 
signals like those due to pendulous motion, magnetometer electronic noise, and A/D converter quantiza- 
tion. In addition, it desensitizes the scheme to bias errors and low frequency errors, like those occurring 
at orbital frequency (0.00018 Hz). These are bound to exist in estimating the Earth's magnetic field. 
Hence, the steady-state outputs of the bandpass filters should better reflect the right hand sides in equa- 
tion (7) than the inputs. Inferences about satellite skiprope motion will be drawn from them. 


Two final comments about implementation are appropriate. First, the smaller the £ ot the 
bandpass filters, the better they attenuate unwanted signals, but the longer it takes their outputs to reach 
steady state. The latter can be improved by choosing the initial states of each filter so that one equals the 
initial input to the filter and the other equals zero. Doing this. 4 = 0.1 was found by simulation to give a 
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good compromise between filtering and speed. The other point to make is that the scheme will normally 
require two iterations, since the skiprope frequency can only be estimated a priori. For the first iteration, 
(o„ should be set to the best guess for w SR . At steady state, the filtered output for B SK in equation (7) will 
reveal u> SR more precisely. Then, u)„ can be tuned to it and the process repeated for the second iteration. 


IV. SIMULATION RESULTS 


A computer simulation was developed to test the scheme described in section III, subject to the 
following conditions. (See the appendix for a listing of the program.) The tethered satellite was assumed 
to be in a 300-km orbit with a 28.5° inclination. Its motion relative to the orbital coordinate frame 
described in section II was given by 


Otv3 = O-SMi 

a S2 “ 0i S2B + GL SR2\f sin ( + CLp2M sin tt/4) 

a VI = a .vi/? + a .s7?iA/ COS + OLp).\f sin (<Mpt + tt/4) , 


where 


«V3B = 5° 


CX S2B 


a .V \B ~ I ° 
° i SR2M = 0 . 5 ° 

a .S7? 1 M = 0.75° 

(X P2M ~ 2 ° 

am/ = 2° 


a>.s 7 ? = 277(0.005) rad/s 
top = 277(0.03125) rad/s . 


Hence, the satellite showed skiprope and pendulous motion, plus bias errors. The bandpass filters were 
set so that 
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0) 


0I.V* = 2ir(0.005Hz) 


a = 0.1 . 


which implies that the litters were already tuned to the skiprope frequency. The Earth's magnetic held 
was modeled by a six-displaced-dipole model. The estimated held was given by 


B = 1.15 

(3/7 = Pay. + 5° 
P.*/ = Pa/ + 5° . 


where (B.$h ,$.\y) are its true values. The magnetometer axes were assumed to be aligned with the satel- 
lite axes with measurements made 16 times per second and quantized to an LSB = 0.02 milli-Gauss |6|. 
This corresponds to a 16-bit A/D converter scaled to a range of ±0.65536 Gauss. 

For these conditions, the simulation results are shown in figures 5. 6. and 7. Using the rules out- 
lined in section II for interpreting the plots in figure 7, the following conclusions are drawn about the 
skiprope motion. In tisure 7. observe that B $ if leads Bsic by tt/2 rad. hence. From B yx - if 

be seen that T SR = 200s and so u> s7? = 2 tt(0.005) rad/s. The amplitude o\ B SK shows that & VM = 0.01 25 
rad = 0.72°. Flence. the deflection of the tether^at its midpoint due to skiprope is d = 0.0125 Lrn by 
virtue of equation ( 1 ). At t = 100.300.... s. the T s , axis lies in the plane of i H and B. tilted toward B. At t 
= 200.400.... s. it lies in the samejtlane. but tilted away from B. The tether is -it rad out of phase with 
the projection ol /yi onto the i( -i\ plane. 

In all aspects, the estimated skiprope motion matches the actual reasonably well, in spite of the 
demanding test conditions. This verifies the scheme and demonstrates its robustness. 


V. FINAL COMMENTS 


This paper has presented a unique scheme for reconstructing tethered satellite skiprope motion by 
bandpass filtering satellite magnetometer outputs telemetered to the ground. It was tested in a computer 
simulation with a demanding set of test conditions. This showed it to be quite robust. 

As a final remark, this scheme is not just limited to the tethered satellite skiprope problem posed 
here. Indeed, it has potential application wherever: 

I A body cones about a known axis while measuring a known vector in body-lixed axes. 

2. There is a need to know the time history, or perhaps just the fundamental characteristics ot. 
the conins motion. 
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ODD 


C THIS PROGRAM IS NAMED TSSMAG .FORT ( CRAY6) . IT TESTS SCHEME TO 
C RECONSTRUCT TETHERED SATELLITE ATTITUDE MOTION DUE TO TETHER 
C SKIPROPE USING MAGNETOMETER OUTPUTS. BANDPASS FILTERS ARE USED TO 
C RECONSTRUCT SKIPROPE FREQUENCY COMPONENTS. LOWPASS FILTER IS USED TO 
C RECONSTRUCT AS3 . FOURTH ORDER RUNGA-KUTTA IS USED TO SOLVE FILTER 
C DIFF. EQS. THIS PROGRAM WAS LAST REVISED ON FEB 19*1991. 

REAL NU2 . L A2 * L A3 * GA2 
DIMENSION RVG<3) * BMAG ( 3^ * AA ( 38 ) 

C SPECIFY AND COMPUTE CONSTANTS 

READ (4*2) IDUM 
2 FORMAT (14) 

PI=3. 1415926 
TWOP 1=2 . 0*PI 
PI 02=P I / 2 . 0 


RFD=0. 01 7453292 
DFR= 1 . Q/RFD 

C ALT =ORB IT ALT IN KM 
C REA=RAD IUS OF EARTH IN KM 
GME = 39860 1 .2 
ALT=300 . 0 
RE A=6371 .2 
RV=REA+ALT 

C OMO=ORB I T RATE IN RAD/SEC 
OMO=SQRT (GME/ ( RV**3) ) 

TO=TWOP I /OMO 

LA3=0RBI T I CL I NAT I ON IN RAD 

OMA=ORBIT REGRESSION RATE IN RAD/SEC (SEE THOMPSON'S 
DYNAMICS* P.99 FOR EQUATIONS GIVEN BELOW) 
LA3=RFD*28 . 5 
X J = 1 .637E-03 

PSI=TWOPI*XJ*COS (LA3 ) * (REA/RV) **2 


"INTRO TO SPACE 


OMA=-PSI /TO 

C OMG=EARTH RATE IN RAD/SEC 

0MG=RFD*360. 0/ (24.0*3600.0) 
ASB3=RFD*5 . 0 
ASB2=RFD* 1 .0 
ASB1 =RFD*1 .0 
W5R=TW0PI*0. 0050 
ASR2M=RFD*0 . 50 
ASR1M=RFD*0. 75 
WP=TW0PI*0. 03125 


AP2M=RFD*2 . 0 
AP 1M=RFD*2 . 0 


PH IP=RFD*45 . 0 
WNB=TWOPI*0.0001 
ZETAB=0 . 707 
WNF=TWOP 1*0. 0050 


2ETAF=0 . 100 


DT= 1 .0/16.0 


NPMAX= 16 
TMAX=6000 . 0 
PRINT* *0M0 *0MA *OMG 

C COMPUTE GAIN AND PHASE OF BANDPASS FILTER TUNED FOR 
C SKIPROPE FREQUENCY 
RATI 0=WSR/WNF 
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no non n noon 


XR= 1 . 0-RAT I 0*RAT 10 
XI=2.0*ZETAF*RATI0 
GM=RATI0/SQRT (XR*XR+XI*XI) 

GP=F'I02-ATAN2(XI ,XR) 

PR INT* , ZETAF , WNF , WSR , GM , GP 

C SPECIFY ERRORS IN ESTIMATING EARTH’S MAGNETIC FIELD 
BHEF=0 . 1 0 
BELHE=RFD*5 . 0 
BAZHE=RFD*5 . 0 

C MAGNETOMETER OUTPUT QUANTIZATION PARAMETERS 
NB ITS= 1 6 
BSMAX=0. 65536 
IQ = 1 

IQFST 1 =1 
IQFST2=1 
IQFST3= 1 

C DUMMY VARIABLE FOR PLOT VARIABLES NOT USED 
DUM=0 . 0 

C SPECIFY INITIAL CONDITIONS 

T=0 . 0 
NP=NPMAX 
START= 1 . 0 

GA2=ANGLE FROM VERNAL EQUINOX TO PRIME MERIDIAN, NORTH IS POSITIVE, 
IN RAD 

LA2=L0NGI TUDE OF ASCENDING NODE IN RAD 

NU2=P0SIT ION OF VEHICLE IN ORBIT WRT ASCENDING NODE IN RAD 
GA2=0 . 0 
LA2=0. 0 
NU2=0 . 0 

INITIAL STATES OF FILTERS: 

BS 1PF=0 . 0 
BS1PFD = 0 .0 
BS2PF=0 . 0 
BS2PFD=Q . 0 
BS3PF=0 . 0 
8S3PFD=0 . 0 

I FLGMF=0 GIVES CONSTANT EARTH MAG FIELD WRT LOCAL VERTICAL AXES 
I FLGMF= 1 GIVES 6-DI SPLACED-DI POLE MODEL OF EARTH MAG FIELD WITH 

REALISTIC TETHERED SATELLITES ORBIT MECHANICS 

IFLGMF=1 
10 CONTINUE 

DETERMINE EARTH’S MAG FIELD FLUX DENSITY VECTOR IN GAUSS IN LOCAL 

VERTICAL AXES 

0E11=C0S (NU2)*C0S(LA3) *COS (LA2 ) -SI N <NU2> *SIN (LA2) 

0E12=C0S (NU2 ) #SI N ( LA3) 

0E13=-C0S (NU2) *C0S (LA3 ) #SI N ( LA2) -SIN (NU2 ) *COS ( LA2) 

0E21=-SIN(LA3) *C0S (LA2) 

0E22=C0S (LA3 ) 

0E23=SIN (LA3)*SIN(LA2) 

0E31=SIN (NU2)*C0S(LA3) *CQS (LAE ) +COS (NU2) *S IN (LA2 ) 

0E32=S IN (NU2)#SIN(LA3) 

0E33=-SI N ( NU2) *COS (L A3 ) #SI N ( LA2) +COS (NU2 ) *COS ( LA2 ) 

SGA2=S I N < -GA2) 

CGA2=C0S (-GA2) 
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AEG 1 1 = CG A2 
AEG1 2=0 . 

AEG1 3=-SGA2 
AEG2 1 =0 . 

AEG22= 1 . 0 
AEG23=0. 

AEG3 1=SGA2 
AEG32=Q . 

AEG33=CGA2 

C SUBROUTINE CALL STATMENTS 
TEMBF=GA2-LA2 
CTEM=C0S (TEMBF) 

CLA3=C0S (-LA3) 

SNU2=SIN (-NU2) 

STEM=S IN (TEMBF ) 

CNU2=C0S (-NU2) 

SLA3=SIN (-LA3) 

RVG( 1) =(-SNU2*CTEM*CLA3-STEM*CNU2) *RV 
RVG ( 2 ) = (SNU2*SLA3)*RV 

RVG (3) = ( -SNU2*STEM*CLA3+CTEM*CNU2) *RV 
CALL BF I ELD ( RVG » ST ART r BMAG ) 

C 

BM AGT= SORT (BMAG ( 1 ) **2+BMAG (2 ) **2+BMAG ( 3) **2> 

BEG1 =AEG 1 1 #BMAG ( 1 ) +AEG 12*BMAG ( 2) +AEG13*BMAG (3) 
BEG2=AEG21*BMAG( 1) +AEG22*BMAG(2) +AEG23*BMAG(3) 
BEG3=AEG31 *BMAG ( 1) +AEG32*BMAG ( 2) +AEG33*BMAG(3) 

C 

BU=0E1 1*BEG1 +0E12*BEG2+0E13*BEG3 
BV=0E2 1*BEG1 +0E22*BEG2+0E23*BEG3 
BW=0E31*BEG1+0E32*BEG2+0E33*BEG3 
C I FLGMG=0 GIVES CONSTANT EARTH MAG FIELD IN LOCAL 
C I FLGMG= 1 GIVES VARIABLE EARTH MAG FIELD IN LOCAL 
IF ( I FLGMF . EQ . 1 ) GO TO 15 
B=Q . 3 

8EL=RFD*30 . 0 
BAZ=RFD*30 .0 
BU=-B*COS (BEL) *SIN (BAZ ) 

BV=B*COS (BEL) #COS(BAZ) 

BW=B*SIN (BEL) 

15 CONTINUE 

B=SQRT (BU**2+BV**2+BW**2) 

BP=SQRT (BU*#2+BV**2> 

BUN=-BU 

BEL=ATAN2(BW,BP) 

BAZ=ATAN2(BUN»BV) 

C DETERMINE EULER ANGLES RELATING SATELLITE AXES TO 

C AXES 

THSR=WSR*T 

THP=WP*T+PHIP 

AS3=ASB3 

ASR2=ASR2M*SIN (THSR) 

ASR 1 =ASR 1M*C0S (THSR) 

AP2 = AP2M»S IN (THP ) 

AP 1 = AP 1M*S IN (THP) 


VERTICAL AXES 
VERTICAL AXES 


LOCAL VERTICAL 
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AS2=ASB2+ASR2+AP2 

AS1=ASB1+ASR1+AP1 

C DETERMINE MAGNETOMETER OUTPUTS IN SATELLITE AXES 
SI =S IN (AS1 ) 

Cl=COS (ASI ) 

S2=S IN ( AS2 ) 

C2=C0S ( AS2 ) 

S3=SIN (AS3 ) 

C3=C0S ( AS3 ) 

T 1 1=C2*C3 
T12=C2*S3 
T13=-S2 

T21=S1 *S2*C3-C1*S3 
T22=S1*S2*S3+C1*C3 
T23=S1*C2 

T31=C1*S2*C3+S1*S3 

T32=C1*S2*S3-S1*C3 

T33=C1*C2 

BS1=T1 1 *BU+T 12*BV+T13*BW 
BS2=T21*BU+T22*BV+T23*BW 
BS3=T31*BU+T32*BV+T33*BW 

C MAGNETOMETER OUTPUT QUANTIZATION 

CALL QUANT (IQ, IQFST1 , NB I TS , BSM AX , BS1 ,BS1Q> 

CALL QUANT (IQ, IQFST2 , NBI TS , BSMAX , BS2 , BS2Q) 

CALL QUANT (IQ, IQFST3 , NBI TS , BSMAX , BS3 , BS3Q) 

BS 1QE=BS 1Q-BS1 
BS2QE=BS2Q-BS2 
BS3QE=BS3Q-BS3 

C COMPUTE INPUTS TO FILTERS 
BH= ( 1 . O+BHEF ) *B 
BELH=BEL+BELHE 
BAZH=BAZ+BAZHE 
SBELH=S I N ( BELH ) 

CBELH=COS ( BELH ) 

SBAZH=SIN(BAZH) 

CBAZH=COS(BAZH) 

BS0F'=(BS1Q*CBAZH + BS2*SBAZH) / (BH*CBELH) 

BS1P= (BS1Q+BH*CBELH*SBAZH) /BH 
BS2P=(BS2Q-BH*CBELH*CBAZH> /BH 
BS3F‘= < BS3Q-BH*SBELH) / (BH*CBELH) 

C SPECIFY BETTER INITIAL CONDITIONS FOR FILTER OUTPUTS AT T=0.0 

IF (T .GT . 0 . 0) GO TO 40 
BSOF'F=BSOP 
BS 1PF=BS IP 
BS2PF=BS2P 
BS3PF=BS3P 
40 CONTINUE 

C TEST FOR TIME TO STOP OR STORE DATA IN ARRAYS FOR PLOTTING LATER 
IF (T.GT.TMAX)GO TO 100 
20 IF (NP.LT.NPMAX)GO TO 30 
BELU=DFR#BEL 
BAZU=DFR*BAZ 
BELHU=DFR*BELH 
BAZHU=DFR*BAZH 
AS3H=BS0PF 
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BS 1 CH= (B31PFD/WNF) /GM 

BS2CH= (BS2PFD/WNF) /GM 

BS3CH= (BS3PFD/WNF) /GM 

AA (01) =T 

AA (02) =B 

AA (03) =BELU 

AA (04) =B AZU 

AA (05) =BH 

AA (06) =BELHU 

AA (07) =BAZHU 

AA (08) =BS 1 

AA (09) =BS2 

AA (10) =ES3 

AA ( 1 1 ) =BS1 Q 

AA ( 12) =BS2Q 

AA ( 13) =BS3Q 

AA (14) =BS1QE 

AA (15) =BS2QE 

AA (16) =BS3QE 

AA (17) =BS0P 

AA (18) = BS 1 P 

AA ( 1 9 ) = B52P 

AA (20) =BS3P 

AA (21 ) =BS0PF 

AA (22) =BS 1 PF 

AA (23) =BS2PF 

AA (24) =BS3PF 

AA (25) =BS0PFD 

AA (26) =BS1 PFD 

AA (27) =BS2PFD 

AA (28) =BS3PFD 

AA (29) =BS 1 CH 

AA (30) =BS2CH 

AA (31) =BS3CH 

AA (32) =AS3H 

AA (33) =A SR 1 

AA (34) =ASR2 

AA (35) =AP 1 

AA (36) =AP2 

AA (37) =AS 1 

AA (38) =AS2 

WR I TE ( 7 ) AA 

NP = 0 

30 CONTINUE 

C UPDATE FILTER OUTPUTS, ORBIT/EARTH ANGLES, TIME, AND COUNTER FOR 

C STORING DATA TO PLOT 

CALL LPF2 ( DT , WNB , ZETAB , BSOP , BSOPF , BSOPFD ) 

CALL LPF2 ( DT , WNF , ZETAF ,BS1 P , BS 1PF , BS1PFD) 

CALL LPF2 ( DT , WNF , ZETAF ,BS2P , BS2PF, BS2PFD) 

CALL LPF2(DT , WNF , ZETAF , BS3P , BS3PF » BS3PFD ) 

GA2=GA2+0MG*DT 

LA2=LA2+0MA*DT 

NU2=NU2+0M0*DT 

IF (NU2 .GT. TWOPI) NU2=NU2-TW0P I 
T=T+DT 
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NP=NP+ 1 
GO TO 10 
100 CONTINUE 
STOP 
END 

SUBROUTINE LPF2 ( DT , WN , ZETA , U , X 1 , X2) 

C THIS SUBROUTINE UPDATES STATES FOR 2ND ORDER LOW PASS FILTER USING 

C 4TH ORDER RUNGA-KUTTA 

X1P=X1 
X2P=X2 
SUM1 =0 . 0 
5UM2=0 .0 
DC 5 1=1 ,4 
WFACT= 1 . 0 
DTP=DT/2.0 

IF ( I .EQ. 2 . OR . I .EQ.3) WFACT=2.0 
IF (I ,EQ. 3) DTP=DT 
X 1 D= X2P 

X2D=WN* ( WN* ( U-X 1 P) -2 .0*ZETA*X2P) 

SUM1 =S'JM 1 + WFACT*X1 D 
SUM2=SUM2+WFACT*X2D 
XI P= X 1 + DTF‘*X ID 
X2P=X2+DTP*X2D 
5 CONTINUE 

Xl=Xl+SUMl*DT/6.0 
X2 = X2+SUM2*DT/6. 0 
RETURN 
END 

SUBROUTINE BF I ELD ( RV , START , BMAG ) 

C THIS IS 6-DISPLACED-DIPOLE MODEL OF THE EARTH'S MAG FIELD 

C FLUX DENSITY IS GIVEN IN GAUSS 

DIMENSION RV(3).R0(12f3) »R(3) * BMAG (3) , EM (12. 3) 

IF (START .NE. 1 . ) GO TO 4 
RO (1 , 1 5 = 500.805 
R0(2,l>=-292.9318 
RO (3,15=455.075 
RO (4,15=414.636 
RO (5 , 1 5 =-161 1.292 
RO (6,15 = 14.958 
C 

RO <1 ,25 = 146.300 
RO (2,2) = -818 .6409 
RO (3 , 2 5 = -902 .802 
R0<4, 25 = 1540. 2948 
R0(5, 25=1095. 094 
R0(6, 25=1243. 2935 
C 

R0(1 ,35=713.994 
R0(2, 35=583. 55 
R0(3, 35=791. 483 
RO (4 , 3 5 =-65 . 87 8 
R0(5,3)=-232.390 
R0(6,3)=-849.084 
C 

EM ( 1 , 1 5 = 1 . 423489E6 
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c 


EM <2 , 1 ) =4 . 558039E6 
EM (3 , 1 ) = ~7 .045997E6 
EM (4 , 1 ) = - . 5556 10E6 
EM (5 > 1 ) = . 1 54356E6 
EM (6 , 1 ) =2 . 9203 14E6 

EM (1 »2)=-13.247848E6 
EM (2,2)=3.575713E6 
EM (3 ,2 ) =3. 744393E6 
EM (4 >2 ) =- . 696502E6 
EM(5,2)=-.593411E6 
EM (6,2)=4.470056E6 

EM ( 1 , 3 ) =2 . 104301 E6 
EM (2,3>=-3 .079872E6 
EM (3,3)=-4. 175109E6 
EM (4 ,3 ) =6. 80 1960E6 
EM (5 i 3 ) =- . 856652E6 
EM (6 ,3) =-l . 303223E6 

ST ART = 0. 0 

4 DO 1 J = l,3 

1 BMAG ( J ) =0 . 

DO 3 1=1 , 6 
DO 2 J = i .3 

2 R(J)=RV(J) -ROU, J) 

RM=R ( 1 )*R< 1) +R (2)*R<2> +R(3)*R(3) 

B=3 . * ( R ( 1) #EM ( I» 1) +R (2 ) *EM ( I »2)+R(3)*EM(I»3> )/RM 
RM=RM**1 .5 
DO 3 J=1 ,3 

3 BMAG ( J) =BMAG ( J) - (EM ( I , J) -B*R ( J) ) /RM 

C CONVERT FLUX DENSITY FROM WEBERS/ ( METERS**2> TO GAUSS 
DO 5 J=1 ,3 

5 BMAG (J)=BMAG ( J > * 1 . 0E+04 
RETURN 

END 

SUBROUTINE QUANT (IQ, IQFST , NB ITS, XMAX , X , XQ) 

IF (IQFST .NE. 1) GO TO 10 
XLSB= (2. 0*XMAX ) / (2 .0**NBITS) 

XMAXP=XMAX-XLSB 
XM I N=- XMAX 

PRINT#, NBITS, XMAX, XL SB, XMAXP,XMIN 
I QFST=0 
10 CONTINUE 
XQ = X 

IF (IQ.EQ.O)GO TO 20 
XDUM = X Q 

IF (XDUM.LT .XMIN) XDUM=XMIN 

IF (XDUM.GT .XMAXP) XDUM=XMAXP 

XDUM=X DUM+ XMAX +XLSB/ 2.0 

XDUM = X DUM/ XLSB 

IXDUM=XDUM 

XDUM=I XDUM 

XDUM = X DUM* XLSB 
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XQ = X DUM- XM AX 
20 CONTINUE 
RETURN 
END 
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